# **************************** LICENSE START ***********************************
#
# Copyright 2012 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************

extern gradientb(f:fieldset) "fortran90" inline
!
!  "GRADIENTB" COMPUTES THE GRADIENT OF A FIELD
!
!  THIS PROGRAM IS A MODIFIED VERSION OF THE FORMER
!  "GRADIENT" TO TAKE INTO ACCOUNT THE UNITS.
! 
!  THE UNITS ARE IN THE INTERNATIONAL SYSTEM
!
!  GRIBEX version:    October,  1996
!  GRIB_API version:  March,    2010
!  MFI version:       November, 2010

  PROGRAM GRADIENTB
  USE grib_api
  IMPLICIT NONE

  INTEGER  fieldset_in, fieldset_out, icnt
  INTEGER  grib_id, isize, istatus, i
  INTEGER  byte_size
  REAL*8,  ALLOCATABLE :: grib_in(:)
  REAL*8,  ALLOCATABLE :: grib_out_u(:)
  REAL*8,  ALLOCATABLE :: grib_out_v(:)

                                    !-- GET FIRST ARGUMENT AS A FIELDSET.
                                    !-- icnt IS THE NUMBER OF FIELDS
  CALL mfi_get_fieldset( fieldset_in, icnt )

                                    !-- CREATE A NEW OUTPUT FIELDSET
  CALL mfi_new_fieldset( fieldset_out )

                                    !-- LOOP ON FIELDS
  DO i=1, icnt
                                    !-- GET NEXT FIELD FROM INPUT FIELDSET
     CALL mfi_load_one_grib( FIELDSET_IN, grib_id )

                                    !-- ALLOCATE ARRAYS, GET FIELD VALUES
     CALL grib_get_size( grib_id, 'values', isize )
     ALLOCATE( grib_in(isize), grib_out_u(isize), grib_out_v(isize) )

     CALL grib_get_real8_array( grib_id, 'values', grib_in, istatus )

                                    !-- VALIDATE AND DERIVE OUTPUT
     CALL valid( grib_id )
     CALL grad( grib_in, grib_out_u, grib_out_v )

                                    !-- SET OUTPUT AS U COMPONENT OF WIND
     CALL grib_set_real8_array( grib_id, 'values', grib_out_u, istatus )
     CALL grib_set_int( grib_id, 'paramId', 131 )

                                    !-- ADD IT TO THE OUTPUT FIELDSET
     CALL mfi_save_grib( fieldset_out, grib_id )

                                    !-- SET OUTPUT AS V COMPONENT OF WIND
     CALL grib_set_real8_array( grib_id, 'values', grib_out_v, istatus )
     CALL grib_set_int( grib_id, 'paramId', 132 )

                                    !-- ADD IT TO THE OUTPUT FIELDSET
     CALL mfi_save_grib( fieldset_out, grib_id )

                                    !-- RELEASE MEMORY
     CALL grib_release( grib_id )
     DEALLOCATE( grib_in, grib_out_u, grib_out_v )

  END DO
                                    !-- RETURN THE RESULT
  CALL mfi_return_fieldset( fieldset_out )

  STOP
END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--
!--  USER ROUTINE TO CHECK VALIDITY OF INPUT FIELD
!--  VALID FOR A GLOBAL FIELD, LAT/LONG, 1.5 DEG GRID
!--

SUBROUTINE valid( grib_id )
  USE grib_api
  INTEGER grib_id
  INTEGER ivalue
  REAL*8  rvalue

  CALL grib_get_int(grib_id, 'dataRepresentationType', ivalue)
  IF( ivalue .NE. 0 ) CALL mfi_fail("GRID not lat/lon")

  CALL grib_get_real8(grib_id, 'iDirectionIncrementInDegrees', rvalue)
  IF( rvalue .NE. 1.5 ) CALL mfi_fail("GRID not 1.5/1.5")

  CALL grib_get_real8(grib_id, 'jDirectionIncrementInDegrees', rvalue)
  IF( rvalue .NE. 1.5 ) CALL mfi_fail("GRID not 1.5/1.5")

  CALL grib_get_int( grib_id, 'Ni', ivalue )
  IF( ivalue .NE. 240 ) CALL mfi_fail("GRID not global")

  CALL grib_get_int( grib_id, 'Nj', ivalue )
  IF( ivalue .NE. 121 ) CALL mfi_fail("GRID not global")

  RETURN
END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--
!--  DERIVE GRADIENT OF INPUT FIELD F (BENITO ELVIRA, IM)
!--  FA = HORIZONTAL GRADIENT, FB = VERTICAL GRADIENT
!--

SUBROUTINE GRAD (F, FA, FB)

  !-- DIMENSIONS CORRESPONDING TO 1.5 x 1.5 GRID
  DIMENSION F(240,121), FA(240,121), FB(240,121)

  PI = ACOS(-1.0)
  RT = 6371000.0
  CB = (RT*1.5*PI)/180.0
                                    !-- COMPUTE HORIZONTAL GRADIENT
  DO I = 1, 121

     C = COS( (90.0-I*1.5 + 1.5)*PI/180.0 )
     FA(1,i) = (F(2,i)-F(240,i)) / (2.0*C*CB)
     FA(240,i) = (F(1,i)-F(239,i)) / (2.0*C*CB)

     DO J = 2, 239
        FA(j,i) = (F(j+1,i)-F(j-1,i)) / (2.0*C*CB)
     END DO

  END DO
                                    !-- COMPUTE VERTICAL GRADIENT
  DO I = 1, 240

     FB(i,1) = 0
     FB(i,121) = 0
     DO J = 2, 120
        FB(i,j) = (F(i,j+1)-F(i,j-1)) / (-2.0*CB)
     END DO

  END DO

  RETURN
END

end inline





# set the area we wish to retrieve data from
#          N,   W,  S,  E

area_xx = [70, -45, 10, 85]


# Retrieve the specific humidity
q = retrieve(
		date	:	-1,
		param	:	"q",
		level	:	700,
		grid	:	[1.5,1.5]
		)

# Get the u and v components of the wind
u = retrieve(
		date	:	-1,
		param	:	"u",
		level	:	700,
		area	:	area_xx,
		grid	:	[1.5,1.5]
		)
v = retrieve(
		date	:	-1,
		param	:	"v",
		level	:	700,
		area	:	area_xx,
		grid	:	[1.5,1.5]
		)


# Compute the gradient of Q
q = gradientb(q)

# Extract the area we are calculating on
q = read ( area : area_xx, data : q)



# Compute the advection of Q
a = q[1]*u + q[2]*v
a = -a * (10 ^ 8) # units will be 10e-8 (kg/kg)/sec


# Plot positive advection in blue, negative in red
contour_common = (
		contour_level_selection_type : "interval",
		contour_interval             : 3,
		contour_label                : "on",
		contour_label_height         : 0.25,	
		contour_highlight            : "off",
		contour_hilo                 : "on",
		contour_hilo_type            : "number",
		contour_hilo_format          : "F5.1",
		contour_hilo_height          : 0.3
		)

cont_n = mcont(
		contour_common,
		contour_max_level       : -0.0001,
		contour_line_colour     : "red",
		contour_label_colour    : "red",
		contour_lo_colour       : "red"
		)

cont_p = mcont(
		contour_common,
		contour_min_level       : 0.0001,
		contour_line_colour     : "blue",
		contour_label_colour    : "blue",
		contour_hi_colour       : "blue"
		)

# A plot window
acoast = mcoast(
		map_coastline_resolution        : "low",
		map_grid_longitude_increment    : 10,
		map_coastline_land_shade        : "on",
		map_coastline_land_shade_colour : "cream"
		)

ps_atlantic = geoview(
		map_projection  : "polar_stereographic",
		area            : [30,-25,50,65],
		coastlines      : acoast
		)

page = plot_page(
		view            : ps_atlantic
		)

dw = plot_superpage(
		custom_width    : 29.7,
		custom_height   : 21,
		pages           : page
		)



# Now plot the result

plot(dw,a,cont_p,cont_n)

